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Gutzwiller functions are popular variational wavefunctions for correlated electrons in Hubbard 
models. Following the variational principle, we are interested in the Gutzwiller parameters that 
minimize e.g. the expectation value of the energy. Rewriting the expectation value as a rational 
function in the Gutzwiller parameters, we find a very efficient way for performing that minimization. 
The method can be used to optimize general Gutzwiller-type wavefunctions both, in variational and 
in fixed-node diffusion Monte Carlo. 
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I. INTRODUCTION 

The Hubbard HamiltoniarJjHl is a basic model for 
studying correlated electrons. Despite its simple form 
it is very difficult to solve, the main complication arising 
from the fact that the interaction term is diagonal in real 
space, while the kinetic energy is simple in momentum 
space. Except for special cases, namely the oneta and the 
infinite-dimensionaH Hubbard model, no exact solutions 
are known. For other dimensions we thus have to use ap- 
proximate methods. The variational method has the ad- 
vantage of being applicable over the whole range of inter- 
action strength. Naturally, the quality of variational cal- 
culations depends critically on the trial function and on 
our ability to optimize its parameters. A good trial func- 
tion has to balance the opposing tendencies of the kinetic- 
and the interaction term. Gutzwiller proposed such 
a wavefunction for the_Hubbard model, the Gutzwiller 
wavefunction (GWF).ElCrQ It introduces correlations into 
a Slater determinant by means of a correlation factor 
that is local in configuration space. Like, the Jastrow fac- 
tor in continuum-space wavefunctionsjj it works by re- 
ducing the weight of configurations with large potential 
energy. Unfortunately, also the GWF is hard to treat ex- 
actly. Again, analytical results have only been obtained 
in one dimensionjttj and in infinite dimensions,^ whepa 
the Gutzwiller approximation (GA) becomes exact BE-il 
For 1 < d < oo we have to resort to nuiptpical meth- 
ods like variational Monte Carlo (VMC)Mka Another 
variational method that has been developed recently is 
fixed-node diffusion Monte Carlo (FNDMC) for Fermions 
on a lattice .E3'E£l Like variational Monte Carlo it uses a 
trial wavefunction, but the results are much closer to the 
exact ground state and depend much less on the trial 
wavefunction. Since these Monte Carlo methods work 
in configuration space it is straightforward to implement 
trial functions with more correlation factors. Such gen- 
eralized Gutzwiller-type wavefunctionsLj that include for 
example correlations between empty and doubly occu- 
pied sites giveuimprovements over the original Gutzwiller 
wavefunction.E3 



Our goal is to optimize the parameters in Gutzwiller- 
type wavefunctions in the framework of quantum Monte 
Carlo calculations. General, methods for achieving this 
are correlated samplingL3E3 or th^ recently developed 
stochastic gradient approximation.^ Exploiting the par- 
ticular form of Gutzwiller wavefunctions, we find a differ- 
ent approach, which is equivalent to correlated sampling: 
We observe that expectation values can be rewritten as 
the quotient of two polynomials (i.e. a rational function) 
in the Gutzwiller parameters. Estimating the coefficients 
of these polynomials in a single Monte Carlo run, we 
can then easily minimize the energy expectation value 
by finding the minimum of the corresponding rational 
function. The implementation of this idea in variational 
Monte Carlo is described in section ||. We show how 
the optimization works in practice and give some appli- 
cations. In Sec. Ill the optimization method is adapted 
to work also in fixed-node diffusion Monte Carlo. This 
allows us to study the effect of changing the fixed-node 
constraint in FNDMC. Applications and results for wave- 



functions with more parameters are given in Sec. IV 



II. VARIATIONAL MONTE CARLO 

To be specific, we consider the Hubbard model 



(1) 



where cj a creates an electron with spin a on site i and 

Tii,cr = c la c ia- The sum i n t ne kinetic term is over 
nearest-neighbor pairs, and the hopping matrix element 
is used as the energy scale, i.e. t = 1. The interac- 
tion term in ([!]) may also be written as U D, where 
D = rii^rii^ is the operator that counts the number 
of doubly occupied sites. The Gutzwiller wavefunction is 
then given by 



|*(.g)) = <? D |$o> 
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FIG. 1. Weight of configurations with given number D of double occupancies for the Gutzwiller wavefunction. (The weight 
w(D) is defined as the sum of *t(7?) 

over all configurations R with D doubly occupied sites.) Reducing the Gutzwiller factor 
g suppresses configurations R for which the interaction energy E[ n t(R) — U D(R) is large — at the expense of increasing the 
kinetic energy. The results shown are for 101 + 101 electrons on a 16 x 16 lattice. 



where < g < 1 is the Gutzwiller parameter and |$o) 
is the ground state wavefunction of the non-interacting 
system (U = 0). The Gutzwiller factor builds correla- 
tions into the Slater determinant by reducing the weight 
of configurations with large interaction energy, i.e. large 
number of doubly occupied sites. This is visualized in 
Fig. | 

In the following we briefly review the formalism of vari- 
ational Monte Carlo for lattice systems and discuss a 
way to improve the efficiency of the simulation, especially 
for systems with strong correlation. Then we introduce 
the method for optimizing the Gutzwiller parameter in a 
VMC calculation. 



A. Modified Metropolis algorithm 

To calculate the energy expectation value for the 
Gutzwiller wavefunction we have to perform a sum over 
all configurations R: 

= ,, T , |, T , v = ^ , Tj9 7 > {•$) 



(*t|*t) 



where we have introduced the local energy for a configu- 
ration R: 



E\ oc (R) — 



(*t\R') (R'\H\R) 
(*t\R) 



(4) 



The prime in the last equation indicates that the sum is 
restricted to configurations R' that are connected with 
R by the Hamiltonian, i.e. configurations that can be 
reached from R by hopping one electron to a nearest- 
neighbor site. We call such configurations nearest- 
neighbors in configuration space. 

Since the number of configurations grows combinatori- 
ally with system-size, the sums in (^|) can in general only 
be performed for very small systems. For larger problems 



we can use Monte Carlo. The idea is to perform a ran- 
dom walk in the space of configurations, with transition 
probabilities p(R — ► R') chosen such that the configu- 
rations i?vMC in the random walk have the probability 
distribution function ^^(R). Then 



B VMC — ; ^Ti 



(•5) 



Rvmc 



where the last equality holds within statistical error-bars. 
The transition probabilities are given by p(R — > R') = 



l/N min[l,*|(i?')/*|(i?)], 



N being the maximum 



number of possible transitions .E3 This choice fulfills de- 
tailed balance 



^ 2 T {R)p{R^ R') 



^ 2 T {R')p(R' 



R). 



(6) 



For ergodicity it is sufficient to consider only transi- 
tions between nearest neighboring configurations. The 
standard prescription is then to propose a transition 
R —> R' with probability l/N and accept it with proba- 
bility mm[l^\{R')/^ 2 T {R)]. This works well for U not 
too large. In strongly correlated systems, however, the 
random walk will stay for long times in configurations 
with a small number of double occupancies D(R), since 
most of the proposed moves will increase D an4 hence 



be rejected with probability w 1 - g D ( R ^- D ( R )JB There 
is, however, a way to integrate the time the walk stays 
in a given configuration out. To see how this works, we 
first observe that for the local energy (||) the ratio of the 
wavefunctions for all transitions induced by the Hamil- 
tonian have to be calculated. This in turn means that 
we also know all transition probabilities p(R — > R'). We 
can therefore eliminate any rejection (i. e. make the ac- 
ceptance ratio equal to one) by proposing moves with 
probabilities 



p(R -> R') 



p{R -» R') 



p(R -> R') 



J- Pstay 



(7) 



Checking detailed balance (g) we find that now we are 
sampling configurations Rvmc from the probability dis- 
tribution function 'Fy(-R) (1 — Pstay(R))- To compensate 
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for this we assign a weight w(R) = 1/(1 — Pstay(-R)) to 
each configuration R. The energy expectation value is 
then, within statistical error-bars, given by 



En 



w(R) 



(8) 



The above method is quite efficient since it ensures that 
in every Monte Carlo step a new configuration is visited. 
In other words, the acceptance ratio is one. Instead of 
staying in a configuration where is large, this con- 
figuration is weighted with the expectation value of the 
number of times the simple Metropolis algorithm would 
stay there. This is particularly convenient for simula- 
tions of systems with strong correlations: Instead of hav- 
ing to do longer and longer runs as U is increased, the 
above method produces, for a fixed number of Monte 
Carlo steps, results with comparable error-bars. 



B. Optimization 

We now turn to the problem of minimizing the energy 
expectation value (^|) as a function of the variational pa- 
rameters in the trial function. To this end we could sim- 
ply perform independent VMC calculations for a set of 
different parameters. It is, however, difficult to compare 
the energies from independent calculations since each 
VMC result comes with its own statistical errors. (This, 
problem can be avoided with correlated sampling Jl3'E2l 
The idea is to use the same random walk in calculating 
the expectation value for different trial functions. This 
reduces the relative errors and hence makes it easier to 
find the minimum. 

Let us assume then that we have generated a random 
walk {-Rvmc} for the trial function vpy. Using the same 
random walk, we can also estimate the energy expecta- 
tion value (|5|) for a different trial function ^t- To do so 
we have to compensate for the fact that the configura- 
tions have the probability distribution tyj, instead of 5»y 
by introducing reweighting factors 



En 



Likewise, (||) is reweighted into 

E fivMC w(R) E ioc (R) * 2 t (R)/*t(R) 



Ej 1 



Y,R VMC ™(R)n(R)/n(R) 



(9) 



(10) 



Also the local energy E\ oc (R) can be rewritten to contain 
the new trial function only in ratios with the old one. For 
Gutzwiller functions this implies a drastic simplification. 
Since they differ only in the Gutzwiller factor, the Slater 
determinants cancel, leaving only powers (g/g) D ^: 



E T {g) 



E RvMC EUR) (g/g) 2D{R) 



(11) 




-279.5 



-280.0 



-280.5 
0.50 



0.55 



0.60 



0.65 



FIG. 2. Correlated sampling for the Gutzwiller parameter 
g. The results shown are for 101 + 101 electrons on a 16 x 16 
lattice and (7 = 4. The full curve shows Ex(g) calculated 
using g — 0.65 as Gutzwiller parameter. The predicted min- 
imum <7mm is indicated by the dotted line. The dashed line 
gives the correlated sampling curve obtained from a calcula- 
tion with (? m in = 0.566. Both calculations find the same min- 
imum. The error-bars show the results of independent VMC 
runs for different Gutzwiller parameters. We can clearly see 
the statistical fluctuations in the results of the separate VMC 
runs, which, of course, are absent in the correlated sampling 
curves. 



and 



E loc (R) = -tJ2'(g/g) 



(D(R')-D(R) 



Vt(R') 
^t(R) 



UD{R). 



(12) 



Since the number of doubly occupied sites D(R) for a 
configuration R is an integer, we can then rearrange the 
sums in ( |ll| ) and (|l2] ) into polynomials in g/g. The en- 
ergy expectation value for any Gutzwiller parameter g 
is then given by a rational function in the variable g/g, 
where the coefficients only depend on the fixed trial func- 
tion \9(g)). 

It is then clear how we proceed to optimize the 
Gutzwiller parameter in variational Monte Carlo. We 
first pick a reasonable g and perform a VMC run for 
\^(g)) during which we also estimate the coefficients 
of the above polynomials. We can then easily calcu- 
late Ex{g) by evaluating the rational function in g/g. 
Since there are typically only of the order of a few tens 
non-vanishing coefficients (cf. the distribution of weights 
shown in Fig. |l|), this is a very efficient process. 

Figure || shows how the method works in practice. Al- 
though we deliberately picked a bad starting point we 
still find the right minimum. The correlated sampling 
curves even seem to coincide within the statistical er- 
rors. This is, however, not true for the whole range of 
Gutzwiller parameters. When g differs too much from g, 
the method breaks down. To understand this we again 
turn to Fig. |l]. We see that most configurations in a ran- 
dom walk generated with, say, g — 0.50 will have about 



3 



20 doubly occupied sites. In the Monte Carlo run we 
therefore sample the coefficients for (g/g) 2x20 best while 
the statistics for much larger or smaller powers is poor. 
But it is exactly these poorly sampled coefficients that 
we need in calculating the energy expectation value for 
trial functions with g much different from g. We can 
thus use the overlap of the wavefunctions (ty (g)\ty (g)) 
as a measure for the reliability of the calculated energy 
Ex(g)- Like the energy expectation value itself, it can be 
recast in the form of a rational function, the coefficients 
of which can be sampled during the VMC run: 



m)Mg)) 



(13) 



To get a feeling for the overlap as a function of g for fixed 
g, we can make a rough estimate using the Gutzwillcr 
approximation. Expanding around g we find that the 
overlap looks like a Gaussian: exp[— M (g — g) 2 /a 2 ], with 
M the number of lattice sites. As expected, for g and g 
fixed, the overlap goes to zero exponentially with system 
size. Co is a function of g and the filling. It generally 
decreases with g. This can be understood by looking at 
Fig. |l| as for small g the weights are peaked more sharply 
than for larger Gutzwillcr parameters. For half filling the 
width of the Gaussian is given by do = V^g 2(1 + g). 

The relation between the overlap and the reliability of 
the expression ( pT[ ) can best be seen by comparing VMC 
calculations to exact energies E(g). An example is shown 
in Fig. |^. There we compare the results of VMC calcula- 
tions for finite Hubbard chains of different size with the 
exact result for the infinite chainBEJ Clearly there are 
systematic errors in the energy coming from finite size 
effects in the VMC calculations. But apart from that we 
find a remarkable agreement between the exact energy 
E(g) and the result of the correlated sampling, even for 
fairly small overlaps. It is also evident from the figure 
that the overlap for given g and g decreases with system 
size, making optimizations more and more difficult, the 
larger the system. 

We finally mention some straightforward modifications 
of the scheme we have described above. There are sit- 
uations where it is more appropriate to minimize the 
variance iia — jthe local energy cr 2 (g) rather than the en- 
ergy E(g)£Q Since the variance can also be rewritten in 
terms of a rational function in g/g, variance optimiza- 
tion can be implemented in much the same way as the 
energy minimization that we have described here. Fur- 
thermore, it is clear that the method is not restricted 
to the plain Gutzwiller wavefunction but can be gener- 
alized to trial functions with more correlation factors of 
the type r c ^ . As long as the correlation function c(R) 
is integer-valued on the space of configurations, expecta- 
tion values for such trial functions can still be rewritten 
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FIG. 3. Comparison of the correlated sampling results for 
finite Hubbard chains to the exact result for the infinite chain 
with U = 2. The energies from the VMC calculations are 
indicated by the error-bars in the upper figure. The curves 
give the corresponding results -Evmc(p) from evaluating ex- 
pression @. Clearly there are finite-size errors, especially 
for the small systems, but apart from that the agreement 
between the correlated sampling and the exact energy E(g) 
for the infinite system is remarkable, even for fairly small 
overlaps of the trial functions, as can be seen from the lower 
panel. Eventually, however, when the overlaps become too 
small (here < 0.4) correlated sampling breaks down. 

as rational functions. The only difference to the simpler 
case described above is that now the rational function is 
multivariate, reflecting the fact that there is more than 
one variational parameter. 



III. FIXED-NODE DIFFUSION MONTE CARLO 

We now turn to the optimization of the trial function 
in fixed-node diffusion Monte Carlo (FNDMC). In this 
method the ground state wavefunction of a given Hamil- 
tonian is projected out from a trial function by repeated 
application of a suitable operator. For Fermion systems 
this approach is plagued by the infamous sign-problem, 
which causes an exponential decay of the signal-to-noise 
ratio during a simulation. One approach to evade this 
problem is to use a trial function to fix the nodes of 
the wavefunctions. This is the fixed- node approximation, 
which gives variational estimates of the ground state en- 
ergy as a function of the trial function. Clearly a change 
in the Jastrow factor, will not change the nodes of the 
trial function. Hence one would expect the FNDMC re- 
sults to be independent of such changes, r^his is actu- 
ally true for systems in continuum-space. Eac3 For lat- 
tice problems, the fixed-node approximation is somewhat 
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more subtle, and, as it turps, opt, results also depend on 
the Gutzwiller parameters .11311a Since fixed-node diffusion 
Monte Carlo, like VMC, is a variational method, it is then 
important to optimize the Gutzwiller factors. This opti- 
mization can again be done using correlated sampling. As 
in the preceeding section, we can rewrite the expression 
for the energy as a function of the Gutzwiller parameter 
in terms of a ratio of polynomials. Now, however, the 
order of the polynomials is not a fixed, small number but 
increases with the number of Monte Carlo steps. It is 
therefore not feasible to sample the coefficients directly, 
as we did above. 

We will first briefly review the aspects of fixed- node dif- 
fusion Monte Carlo that are important for implementing 
an optimization scheme for Gutzwiller parameters, then 
describe how to perform the optimization, and finally 
show how the method works in practice. For more in- 
depth discussions and some applications of the FNDMC 
method see, e.g., Refs. 0,p]M0 



A. Fixed-node approximation 

Diffusion Monte CarloEl allows us, in principle, to sam- 
ple the true ground state of a Hamiltonian H . The basic 
idea is to use a projection operator which has the lowest 
eigenstate as a fixed point. For a lattice problem, where 
the spectrum is bounded E n £ [Eo, E max ], the projection 
is given by 



= [1-t(H-E )] |* (n) ), 



(14) 



starting with |*(°)) = |* T ). If r < 2/(E max - E ) and 
|^t) has a non- vanishing overlap with the ground state, 
the above iteration converges to |^o)- There is no time- 
step error involved. Because of the prohibitively large 
dimension of the many-body Hilbert space, the matrix- 
vector product in (|l4|) cannot be done exactly. Instead, 
we rewrite the equation in configuration space 



J2\R'){R'\y {n+1) ) = 

J2\r')(r'\i-t(h 



(15) 



E )\R)(R\^) 



R,R' 



=F{R.',R.) 



and perform the propagation in a stochastic sense: tpw 
is represented by an ensemble of configurations R with 
weights w(R). The transition matrix element F(R',R) 
is rewritten as a transition probability p(R — > R') times 
a normalization factor m(R',R). The iteration ( |l5| ) is 
then performed stochastically as follows: For each R we 
pick a new configuration R' with probability p(R — > R') 
and multiply its weight by m(R', R). Then the new en- 
semble of configurations R 1 with their respective weights 
represents Importance sampling decisively im- 

proves the efficiency of this process by replacing F(R', R) 
with G(R',R) = (V T \R') F(R',R)/{R\^ t ) so that tran- 
sitions from configurations where the trial function is 



small to configurations with large trial function are en- 
hanced. (051) then takes the form 



^|i?')(*T|.R')(i?'l* ( " +1) ) = 

\R') G(R', R) (* T \R) (R\y (n) ) 

R,R' 

i.e. the ensemble of configurations now represents the 
product v{f T \I/< n \ and the transition probabilities are 
given by p(R -> R') = \G(R', R)\/m(R', R) with 
m(R',R) = sign[G(R',R)]J2w'\G(R",R)\ absorbing 
the sign of G(R' , R) and ensuring normalization. After 
a large number n of iterations the ground state energy is 
then given by the mixed estimator 



w _ (* r |ff|*(»)) 



(* T |*(™)) 



Z B E loc (R)w^(R) 



(16) 



with w^(Rn) = n™=i m {Ri> Ri-i)- As long as the evo- 
lution operator has only non-negative matrix elements 
G(R',R), all weights w(R) will be positive. If, however, 
G has negative matrix elements there will be both con- 
figurations with positive and negative weight. Their con- 
tributions to the estimator (|16|) tend to cancel so that 
eventually the statistical error dominates, rendering the 
simulation useless. This is the infamous sign problem. A 
straightforward way to get rid of the sign problem is to 
remove the offending matrix elements from the Hamilto- 
nian, thus defining a new Hamiltonian H g by 



(R'\H cff \R) 





(R'\H\R) 



if G(R',R) < 
else 



(17) 



For each off-diagonal element (R'\H\R) that has been 
removed, a term is added to the diagonal: 

(R\H cB \R) = (R\H\R) + ^t(R'){R'\H\R)/^ t (R). 

R' sf 

This is the fixed-node approximation for lattice Hamilto- 
nians introduced in Ref. [17]. H c s is by construction free 
of the sign problem and variational, i.e. Eq S > Eq, and 
if if ^ t (R')/^t(R) = #o (•#)/*<) (-R) for all R, R' with 
G(R', R) < the method is exact. 

Fixed-node diffusion Monte Carlo for a lattice Hamil- 
tonian thus means that we choose a trial function from 
which we construct an effective Hamiltonian, the ground 
state of which is determined by diffusion Monte Carlo. 
From the equations defining H c g we can then understand 
how the fixed-node results depend on the trial function. 
Clearly, the off diagonal elements ( |l7| ) only depend on 
the sign of the trial function. Since the Gutzwiller term 
is just a non-negative prefactor, a change of g will not 
affect these matrix elements. On the other hand, the di- 
agonal elements (R\H c g\R) can contain ratios of the trial 
function on neighboring configurations. In many cases 
these configurations will differ in the number of doubly 
occupied sites. Then the Gutzwiller terms will not can- 
cel and the diagonal element of the effective Hamiltonian 
will depend on the Gutzwiller factor. 
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B. Correlated sampling 

We are now in the position to describe how we can 
optimize VPy or, equivalently, optimize H c f[ using corre- 
lated sampling. The idea is again to calculate the energy 
for a modified Hamiltonian H e g using a random walk 
generated for the original Hamiltonian H c g. To find the 
reweighting factors involved in the calculation we rewrite 
the mixed estimator 



E, 



(n) 



(* T \Hcs\& n) ) 



(18) 



At first glance it looks like the local energy Ei oc (R) has 
to be calculated for H c s, but keeping track of the sign- 
flip terms in the diagonal of the effective Hamiltonian we 
find 



(19) 



which is just the local energy for the original Hamilto- 
nian. Hence, the local energy can be dealt with as in 
variational Monte Carlo, cf. eqn. (|12|). 

Correlated sampling now means that in the ex- 
pression ( |l8| ) for the mixed estimator we do Monte 
Carlo for the terms n™=i G(Ri, Ri-i) ^^(Rq) instead of 
n£=i G(Ri, R4-1) ^y(ilo)- The reweighting factors are 
thus given by 



n 



gCgi L gi-i) g|(go) 

\G(Ri,Ri-i) *URo)' 



(20) 



Reintroducing the plain (not importance-sampled) 
projection operator F = 1—t(H—Eq), we can rewrite the 
reweighting factor for a simple Gutzwiller wavefunction 
into 



n 



F(Ri, Ri-i) 



- N D(R )+D(R n ) 



\ F{Rii Ri-i) 

with the ratio F(R\ R)/F(R', R) = 1 for R' ^ R, and 



(21) 



F(R, R) 
F(R,R) 



1 -T 


i: f ^{R'\H\R) (}) 


A(R'.R) 

-E 


1 -T 


E ^^{R'\H\R)-E 

_R> sf 





(22) 

with A(R', R) = D(R') - D(R). This expression is again 
(as in variational Monte Carlo) a low order polynomial 
in g/g. The reweighting factor is, however, a product of 
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FIG. 4. Correlated sampling for the Gutzwiller parameter 
g in fixed-node diffusion Monte Carlo. The results are for the 
same system as in Fig. I (101 + 101 electrons on a 16 x 16 
lattice, U = 4). The error-bars show the results of indepen- 
dent FNDMC runs for different Gutzwiller parameters. The 
lines give the results of correlated sampling for trial functions 
with g = 0.44, 0.48, 0.52, 0.56, 0.60, 0.64, and 0.68. Note that 
in fixed-node diffusion Monte Carlo the energy depends much 
less on the trial function than in variational Monte Carlo: 
The energy scale is reduced by a factor 5 compared to Fig. ^, 
although the range of g- values is expanded by a factor of two. 

many such polynomials, therefore the order of the poly- 
nomial representing ( p0| ) increases with the number n of 
Monte Carlo steps. It is therefore not practical to di- 
rectly estimate the ever increasing number of coefficients 
for the reweighting factor. But since we still can easily 
calculate the coefficients for the factors ( p2[ ) we may use 
them to evaluate the mixed estimator E^ (g) in each it- 
eration on a set of predefined values gi of the Gutzwiller 
parameter. 

To see how the method works in practice we look at 
the same system as in Fig. g: 101 + 101 electrons on a 
16 x 16 lattice, with U = 4. Figure |J shows the result 
of independent FNDMC runs and the corresponding cor- 
related sampling results. We see that our method gives 
good predictions of the optimum Gutzwiller parameter, 
even for runs with trial functions that are quite far from 
optimum. Comparing with variational Monte Carlo, we 
first notice that the variational energy in FNDMC is sub- 
stantially lower (by w 2.5 eV). We furthermore see that 
in the fixed-node method there is a pronounced minimum 
in the energy as a function of the Gutzwiller parameter, 
although the energy dependence is much smaller than in 
variational Monte Carlo. This is to be expected since 
in FNDMC only certain features of the trial function 
(namely the ratio of the trial function across the nodes) 
enter the calculation. The situation is, however, different 
from the fixed-node method for wavefunctions defined in 
continuum-space, for which changing the Jastrow factor 
does not change the position of the nodal surfaces and 
hence also does not systematically change the result of 
the fixed-node calculation. We finally notice that the 
optimum Gutzwiller parameter for FNDMC is slightly 
smaller than that obtained from VMC. 
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IV. APPLICATIONS 



We have made practical use of the method for op- 
timizing Gutzwillcr parameters in our investigations of 
the doped Fullerides. In this context we work with a 
Hubbard-like Hamiltonian that describes the conduction 
electrons in the three-fold degenerate t\ u band: 



H 



tinjn' c \na C jn'<j + ^ 



(ij) nn'a 



i (na)<(n' a') 



(23) 



where im.jn' is the hopping integral between orbital n of 
the molecule on site i and orbital n! of molecule j, and 
U is the Coulomb interaction for electrons on the same 
molecule. For more details on the underlying physics 
see Refs. p9|-[33"[ In the Monte Carlo calculations for 
the above Hamiltonian, we use trial functions of the 
Gutzwillcr- type 



\*T(U ,g))=g D MUo)), 



(24) 



where besides the Gutzwiller parameter g we also use 
different types of Slater determinants $. 



A. More Gutzwiller parameters 



sreening for the t± u elec- 
i2l we determine the 



To study the static dielectric i 
trons in the doped Fullerides we determine tne re- 
sponse of the charge density to the introduction of a 
test charge q placed on molecule i q . To describe the test 
charge the term 



H 1 (q) = qUj2> 



(25) 



-VMC 



236.94 
236.92 
236.90 
236.88 
236.86 
236.84 




0.82 



0.28 u 



-FN-DMC 



236.44 
236.43 
236.42 
236.41 




0.82 



0.78 



0.28" 



FIG. 5. Correlated sampling for the parameters g and h 
in the generalized Gutzwiller wavefunction ^t) = g D h N " |<£>), 
cf. eqn. ( p^ ) in variational (upper plot) and fixed-node diffu- 
sion Monte Carlo (lower plot). The plots show the energy 
as a function of the Gutzwiller parameters g and h, both 
as surfaces and contours. The calculations were done for an 
fee cluster of 64 molecules with 96 + 96 electrons (half-filled 
tiu-band), an on-site Hubbard interaction U = 1.25 eV, and 
a test charge of q = 1/4 (in units of the electron charge). 



is added to the Hamiltonian (|23|). In the spirit of 
the Gutzwiller Ansatz we correspondingly add a second 
Gutzwiller factor to the wavefunction (]24|) that reflects 
the additional interaction term qU Ni : 



|*t(s,/0) =g D h N ^\§). 



(26) 



Finding the best Gutzwiller parameters is now a two di- 
mensional optimization problem. Dealing with polyno- 
mials in the two variables g and h, the method of cor- 
related sampling works as straightforwardly as described 
above for the case of a plain Gutzwiller wavefunction. As 
an example, Fig. |{] shows the result of the optimization, 
both in variational and in fixed-node diffusion Monte 
Carlo, for a cluster of 64 Ceo molecules in an fee arrange- 
ment (periodic boundary conditions) resembling K3C60 
with a test charge q = 1/4. In practice we first optimize 
the parameters in variational Monte Carlo. We then use 
the optimum VMC parameters as starting points for the 
optimization in the more time consuming fixed-node dif- 
fusion Monte Carlo calculations. 



B. Variation of Slater determinant 

In the traditional Gutzwillcr Ansatz, the Slater deter- 
minant $ is the ground state wavefunction of the non- 
interacting Hamiltonian. This is, however, not neces- 
sarily the best choice. An alternative would be to use 
the Slater determinant $({/) from solving the interact- 
ing problem in Hartree-Fock approximation. We can even 
interpolate between the two extremes by doing a Hartree- 
Fock calculation with a fictitious Hubbard interaction Uq 
to obtain Slater determinants Q(Uo). Yet another fam- 
ily of Slater determinants <i>(i? s tag) can be obtained from 
solving the non-interacting Hamiltonian with an added 
staggered magnetic field, which lets us control the anti- 
ferromagnetic character of the trial function. Although 
optimizing parameters in the Slater determinant cannot 
be done with the method described in the preceeding sec- 
tions, the cheap optimization of the Gutzwiller factors 
makes it possible to optimize the overall trial function 
without too much effort. 
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1. Staggered magnetic field 

Introducing a staggered magnetic field we can con- 
struct Slater determinants by solving the non-interacting 
Hamiltonian with an added Zeeman term. To be specific, 
we consider K3C60 which has a half-filled ii u -band. Since 
K3C60 crystallizes in an fee lattice, antiferromagnetism 
is frustrated and the definition of a staggered magnetic 
field is not unique. We split the fee lattice into two sub- 
lattices A and B such that the frustration is minimized. 
The Zeeman term is then given by 

H m = H stag ^ sig n M [«*T _ n ii\ ( 27 ) 

i 

with sign(i) = +1 if i G A and — 1 if i G B, It effec- 
tively introduces an on-site energy which, on the same 
site, has opposite sign for the two spin orientations, and, 
for the same spin orientation, has opposite sign on the 
two sublattices. Therefore, hopping to neighboring sites 
on different sublattices involves an energy cost of twice 
the Zeeman energy. The staggered magnetic field thus 
not only induces antiferromagnetic order in the Slater de- 
terminant but also serves to localize the electrons. This is 
reflected in the fact that the optimum Gutzwiller param- 
eter is much larger for Slater determinants constructed 
from a Hamiltonian with large H stag than for paramag- 
netic Slater determinants. Varying H stag then interpo- 
lates between paramagnetic/itinerant and antiferromag- 
netic/localized wavefunctions. 

The energy expectation values for such trial functions 
as calculated in variational Monte Carlo are shown in Fig. 
^. It shows -EVmc as a function of the antiferromagnetic 
correlation 

(s iSi+1 ) = — ^(n 4T - riii) (n n - n n ), (28) 

where the sum is over the N nearest neighbors. (siSj+i) is 
a monotonous function of H stag . For each different value 
of the Hubbard interaction U we find a curve with two 
minima. One minimum is realized for the non-magnetic 
{Hstag — 0) trial function. The energy as a function of U 
scales roughly like E paia oc — (1 — U/U c ) 2 , as predicted by 
the Gutzwiller approximation. The second minimum is in 
the antiferromagnetic/localized region and scales roughly 
like Eaf oc —t 2 /U, as expected. For small U the non- 
magnetic state is more favorable, while for large U the 
localized Slater determinant gives lower variational ener- 
gies. The crossover is at U c ~ 1.50 eV (dotted line) and 
resembles a first order phase transition. 

2. Hartree-Fock 

An alternative method for constructing Slater deter- 
minants is to use the interacting Hamiltonian with the 




FIG. 6. Variational energy -Evmc for trial functions with 
different character. Plotted are the energies (error-bars, lines 
are to guide the eye) for a Hamiltonian describing K3C60 (pe- 
riodic fee cluster of 32 molecules) with Hubbard interaction 
U — 1.0, 1.1 . . . 1.9, 2.0 eV. Instead of the total energies Etot, 
we plot the difference of Etot and the energy in the atomic 
limit (each site occupied by three electrons), so that the re- 
sults for different U can be readily compared. The trial func- 
tions are of the Gutzwiller type. The Slater determinants were 
determined from diagonalizing the non-interacting Hamilto- 
nian (i.e. setting (7 = 0) with a staggered magnetic field 
-ffstag- This field gives rise to an antiferromagnetic correla- 
tion of neighboring spins, which is plotted on the abscissa. 
For U = 1.5 eV (dotted curve) the minima in in the paramag- 
netic and the antiferromagnetic region have about the same 
energy. 



physical Hubbard interaction U replaced by a param- 
eter Uq and solve it in Hartree-Fock approximation. In 
practice this is done by simply doing an unrestricted self- 
consistent calculation for the finite, periodic clusters un- 
der consideration, starting from some charge- and spin- 
density that breaks the symmetry of the Hamiltonian. 
It is well known that Hartree-Fock favors the antiferro- 
magnetic Mott insulator, predicting a Mott transition at 
much too small values U^ F . It is therefore not surpris- 
ing that good trial functions are obtained for values of 
Uo considerably smaller than U. For Uo close to zero 
the Slater determinant has metallic character, while for 
somewhat larger Uq there is a metal-insulator transition. 
Figure [?] shows the energy as a function of Uq for the 
model of K 3 C6o- We find that the results of variational 
Monte Carlo depend quite strongly on the parameter Uo- 
As expected, for fixed Hubbard interaction U there is 
a transition from the paramagnetic region for small Uq 
to a region where the trial function is antiferromagnetic. 
In fixed-node diffusion Monte Carlo energies are overall 
lowered and the dependence on the trial function is much 
weaker. It seems that here mainly the character (param- 
agnetic/antiferromagnetic) of the trial function matters. 
For small U trial functions with small Uo give lower en- 
ergy, while for large U trial functions with larger Uo are 
favorable. The crossover coincides with the Mott trans*-, 
tion, which takes place between U = 1.50 and 1.75 eVJ9 
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FIG. 7. Dependence of variational (VMC) and fixed-node 
diffusion Monte Carlo (FN-DMC) on the trial function. Uo is 
the Hubbard interaction that was used for the Slater determi- 
nant in the Gutzwiller wavefunction *t(#) = &(Uo). 
The results shown here are the energies (relative to the atomic 
limit) for a Hamiltonian that describes K3C60 (32 molecules), 
with U being varied from 1.25 (lowest curve) to 2.00 eV (high- 
est curve). 



V. SUMMARY 

We have presented a convenient and reliable method 
for optimizing Gutzwiller-type trial functions in Monte 
Carlo calculations. The method is based on the obser- 
vation that the expressions for correlated sampling can 
be rewritten in terms of polynomials in the Gutzwiller 
parameters. The method can be used both in variational 
and fixed-node diffusion Monte Carlo. Because of its reli- 
ability and speed it makes optimizing Gutzwiller param- 
eters essentially an automatic process. This is especially 
convenient when dealing with trial functions that have 
several Gutzwiller parameters as in the example on the 
dielectric screening in K^Cqq- Given that optimizing the 
Gutzwiller parameters is quick and easy we can then fo- 
cus on optimizing the Slater determinant. 
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